Setup

Set up workspace, load data, and load required packages.

rm(list=ls(all=TRUE))

results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na.strings="NA")

library("lme4") #linear mixed modeling
library("lmtest") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("car") #levenes tests
library("emmeans") #post-hoc tests
library("lsmeans") #post-hoc tests
library("plotrix") #calculate std.error
library("plyr")  #splitting, applying, and combining data
library("dplyr") 
library("cowplot") #grid plotting and arrange plots
library("effects") #plot effects of modeling
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots, 
library("tidyverse")
library("stats")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("tidyverse")

Approach:

1) Build and run model  
2) Check for normality of residuals  
3) Check for homogeniety of variance of residuals  
4) Look at model summary  
5) Run anova to check for significance of main effects

Analysis:

  1. Environmental Conditions and Manipulations:
    1. Calculate means
    2. Check environmental manipulations with ANOVAs for treatments
    3. Figure 1: Temp over time
    4. Figure 2: pH over time
  2. Urchin Diameters:
    1. Calculate diameter means for:
      1. Day -24 and
      2. Day 1 and check for differences in urchin diameters from pH or
  3. Response Analyses:
    1. Growth (%)
      1. Growth Analysis
      2. Figure 3: Growth Interaction Plot
    2. Calcification Ratio
      1. Calcification Analysis at the Tips
      2. Calcification Analysis at the Bases
      3. Figure 4: Calcification Interaction Plot
    3. Spine Loss (count)
      1. Spine Loss Analysis
      2. Figure 5: Spine Loss Interaction Plot

1. Environmental Conditions and Manipulations:

a. Calculate environmental means

EnviroMean <- 
  results %>% 
    select("Day", "Temperature", "pH", "Temp","pHspec") %>%
    filter(Day >= "1") %>% 
    drop_na(Day) %>% 
    group_by(Temperature, pH) %>% 
    summarise(meanT = mean(Temp, na.rm = T)); EnviroMean
## # A tibble: 4 x 3
## # Groups:   Temperature [2]
##   Temperature pH        meanT
##   <fct>       <fct>     <dbl>
## 1 ambient     acidified  25.0
## 2 ambient     ambient    25.0
## 3 heated      acidified  27.1
## 4 heated      ambient    27.0

b. Check environmental manipulations with ANOVAs for treatments

Environmental <- 
  results %>% 
    select("Day", "Temp", "pH", "Temperature","pHspec","Header", "TankID") %>%
    filter(Day >= "1", Day != "126") #use only days from start of experiment when desired conditions were reached
  

#temp between high and ambient temp
summary(aov(Temp ~ Temperature, data=Environmental)) #temperature was different between high and ambient treatment
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Temperature   1  782.3   782.3   307.6 <2e-16 ***
## Residuals   771 1960.6     2.5                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 355 observations deleted due to missingness
#pH between high and ambient co2
summary(aov(pHspec ~ pH, data=Environmental)) #pH was different between high and ambient co2 treatments
##              Df Sum Sq Mean Sq F value Pr(>F)    
## pH            1 14.012  14.012    5247 <2e-16 ***
## Residuals   627  1.674   0.003                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 499 observations deleted due to missingness
##HEADERS
#temp between headers in high
aov1<- Environmental %>%
  filter(Temperature=="heated")

summary(aov(Temp ~ as.factor(Header), data=aov1)) #not different between headers
##                    Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header)   3    2.0   0.674   0.277  0.842
## Residuals         381  925.7   2.430               
## 179 observations deleted due to missingness
#temp between headers in ambient temp
aov2<- Environmental %>%
  filter(Temperature=="ambient")

summary(aov(Temp ~ as.factor(Header), data=aov2)) #not different between headers
##                    Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header)   3    0.1  0.0318   0.012  0.998
## Residuals         384 1032.8  2.6897               
## 176 observations deleted due to missingness
#pH between headers in acidified
aov3<- Environmental %>%
  filter(pH=="acidified")

summary(aov(pHspec ~ as.factor(Header), data=aov3)) #not different between headers
##                    Df Sum Sq  Mean Sq F value Pr(>F)
## as.factor(Header)   3 0.0243 0.008111   1.785   0.15
## Residuals         316 1.4360 0.004544               
## 248 observations deleted due to missingness
#pH between headers in ambient pH
aov4<- Environmental %>%
  filter(pH=="ambient")

summary(aov(pHspec ~ as.factor(Header), data=aov4)) #not different between headers
##                    Df  Sum Sq   Mean Sq F value Pr(>F)
## as.factor(Header)   3 0.00394 0.0013127   1.906  0.129
## Residuals         305 0.21001 0.0006885               
## 251 observations deleted due to missingness
##TANKS
#temp between tanks in high temp
aov1<- Environmental %>%
  filter(Temperature=="heated")

summary(aov(Temp ~ as.factor(TankID), data=aov1)) #not different between heated tanks
##                    Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID)  11    4.5  0.4077   0.165  0.999
## Residuals         273  675.9  2.4756               
## 279 observations deleted due to missingness
#temp between tanks in ambient temp
aov2<- Environmental %>%
  filter(Temperature=="ambient")

summary(aov(Temp ~ as.factor(TankID), data=aov2)) #not different between ambient tanks
##                    Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID)  11    0.8  0.0687   0.025      1
## Residuals         276  759.8  2.7531               
## 276 observations deleted due to missingness
#pH between tanks in acidified
aov3<- Environmental %>%
  filter(pH=="acidified")

summary(aov(pHspec ~ as.factor(TankID), data=aov3)) #not different between tanks in acidified
##                    Df Sum Sq  Mean Sq F value Pr(>F)
## as.factor(TankID)  11 0.0252 0.002288   0.521  0.888
## Residuals         216 0.9495 0.004396               
## 340 observations deleted due to missingness
#pH between headers in ambient pH
aov4<- Environmental %>%
  filter(pH=="ambient")

summary(aov(pHspec ~ as.factor(TankID), data=aov4)) #not different between headers
##                    Df  Sum Sq   Mean Sq F value Pr(>F)
## as.factor(TankID)  11 0.00516 0.0004686   0.765  0.675
## Residuals         213 0.13045 0.0006124               
## 335 observations deleted due to missingness

c. Figure 1: Temperature over time

tempsummary<-results %>% 
  select("Day", "Temperature", "Temp") %>%
  drop_na(Temp) %>%
  group_by(Day, Temperature) %>%
  mutate(meanT = mean(Temp)) %>%
  group_by(Temperature) %>% 
  mutate(sd = sd(Temp)) %>%
  mutate(se = sd/sqrt(12)) 

tempplot<-ggplot(data=tempsummary, aes(Day, meanT, color = Temperature)) +
    geom_point(size = 2.5, show.legend=FALSE) +
    geom_errorbar(aes(ymin = meanT-se, ymax = meanT+se)) +
    geom_line(size=1.2) +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="Temperature (°C)", breaks = seq(20,30, by = 1)) +
    scale_shape_discrete(name = NULL,
                         labels = c("Ambient", "High"))+
    scale_color_manual(name = NULL,
                       values = c("gray40", "firebrick3"),
                       labels = c("Ambient", "High")) +
    ggtitle("")+
    theme_classic() +
    theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
    theme(legend.margin = margin(0), 
          legend.position = c(.98,.8), 
          legend.justification = c("right", "top"),
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14, face="bold"));tempplot

#ggsave(filename="Figures/20200611/TempFig1.png", plot=tempplot, dpi=300, width=7, height=5, units="in")

d. Figure 2: pH over time

pHsummary<-results %>% 
  select("Day","pH", "pHspec") %>%
  drop_na(pHspec) %>%
  group_by(Day, pH) %>%
  mutate(meanpH = mean(pHspec)) %>%
  group_by(pH, meanpH) %>% 
  mutate(sd = sd(pHspec)) %>%
  mutate(se = sd/sqrt(12)) 
  
  
pHplot<-ggplot(data=pHsummary, aes(Day, meanpH, color = pH)) +
    geom_point(size = 2.5, show.legend=FALSE) +
    geom_errorbar(aes(ymin = meanpH-sd, ymax = meanpH+sd)) +
    geom_line(size=1.2) +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="pH (Total Scale)", breaks = seq(7.5,8.1, by = .1)) +
    scale_shape_discrete(name = NULL,
                         labels = c("Low pH", "High pH"))+
    scale_color_manual(name = NULL,
                       values = c("skyblue3", "gray40"),
                       labels = c("Acidified", "Ambient"),
                       guide = guide_legend(reverse = TRUE)) +
    ggtitle("")+
    theme_classic() +
    theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
    theme(legend.margin = margin(0), 
          legend.position = c(.98,.8), 
          legend.justification = c("right", "top"),
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14, face="bold"));pHplot

#ggsave(filename="Figures/20200611/pHFig1.png", plot=pHplot, dpi=300, width=7, height=5, units="in")

2. Urchin Diameters:

a. Calculate diameter means

i. Day -24

Assemble data for diameters of initial collection of urchins from hatchery

##Initial Collection = Day -24
InitialDiameter <-
  #seperate out the initial sizes of urchins on day -24 after initial collection.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-24") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- std.error(InitialDiameter$Diameter)

InitialMean
## [1] 7.529348
InitialSE
## [1] 0.1068257

Results of linear mixed model effect of diameter (day -24)

InitialMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=InitialDiameter)
anova(InitialMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## Temperature    0.0005035 0.0005035     1    19  0.1144 0.7389
## pH             0.0102269 0.0102269     1    19  2.3232 0.1439
## Temperature:pH 0.0003840 0.0003840     1    19  0.0872 0.7709

Check assumptions of diameter model (day -24)

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(InitialMod)) 

## [1] 27  4
shapiro.test(residuals(InitialMod)) #no pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(InitialMod)
## W = 0.52108, p-value = 4.878e-11
hist(residuals(InitialMod)) #eh

# 2. Equal variances
leveneTest(residuals(InitialMod)~InitialDiameter$Treatment) #No pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.9592 0.04307 *
##       42                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(InitialMod),resid(InitialMod,type="pearson"),col="blue")

ii. Day 1

Assemble Data for diameters on day 1 of experiment when desired conditions were reached

##After ramp up and desired conditions are reached to begin experiment = Day 1

Day1Diameter <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- std.error(Day1Diameter$Diameter)

Results of linear mixed model effect on diameter (day 1)

Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.04177 0.04177     1    19  0.1922 0.66606  
## pH             0.38094 0.38094     1    19  1.7525 0.20128  
## Temperature:pH 0.89969 0.89969     1    19  4.1389 0.05611 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of diameter model (Day 1)

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(Day1Mod)) #not super normal but... see below

## [1]  1 24
shapiro.test(residuals(Day1Mod)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(Day1Mod)
## W = 0.99584, p-value = 0.9998
hist(residuals(Day1Mod)) #looks normalish

# 2. Equal variances
leveneTest(residuals(Day1Mod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.7309 0.5394
##       42
plot(fitted(Day1Mod),resid(Day1Mod,type="pearson"),col="blue")

3. Response Analysis

a. Growth

i. Growth Analysis

Assemble data and calculate means using day 1 as initial diameter and calculate means

## [1] 16.03
## [1] 0.2550688
## # A tibble: 4 x 4
## # Groups:   Temperature [2]
##   Temperature pH         mean  s.e.
##   <fct>       <fct>     <dbl> <dbl>
## 1 ambient     acidified  327. 16.3 
## 2 ambient     ambient    323. 11.0 
## 3 heated      acidified  352. 15.0 
## 4 heated      ambient    376.  7.15

Growth linear mixed model results

#LMM for growth:
modGrow <-
  resultsGrow %>% 
  lmer(Growth~ Temperature * pH + (1|TankID), data=.)

#Need to decide on model type that we want to use. 

summary(modGrow)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 417.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90063 -0.32406  0.08141  0.32903  1.94882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 1929.2   43.92   
##  Residual              283.1   16.83   
## Number of obs: 46, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                  326.595     18.578  19.000  17.580 3.28e-13
## Temperatureheated             25.439     26.273  19.000   0.968    0.345
## pHambient                     -3.424     26.273  19.000  -0.130    0.898
## Temperatureheated:pHambient   27.731     38.073  19.000   0.728    0.475
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.488 -0.690 -0.690
anova(modGrow, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    1169.41 1169.41     1    19  4.1304 0.05634 .
## pH               74.92   74.92     1    19  0.2646 0.61290  
## Temperature:pH  150.20  150.20     1    19  0.5305 0.47527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for growth model

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow)) 

## [1] 43 20
shapiro.test(residuals(modGrow)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow)
## W = 0.95491, p-value = 0.07263
hist(residuals(modGrow)) #looks normal

# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  2.1906 0.1033
##       42
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")

ii. Growth Interaction Plot

b. Calcification Ratio

i. Calcification Ratio Analysis at the tips

Assemble data for chosen spine tips (distal end)

### Model 1: calcification at the tips of spines

resultsTip <- 
  #create data using only the SEM images of the tips of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Tip")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results for spine tips with linear mixed model

#LMM for calcification at the tips of spines 
modTip <-
  resultsTip %>% 
  lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)
## boundary (singular) fit: see ?isSingular
summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 64.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76413 -0.74507 -0.02052  0.63045  2.07630 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.0000   0.0000  
##  Residual             0.1356   0.3682  
## Number of obs: 67, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                  1.69118    0.08930 63.00000  18.938   <2e-16
## Temperatureheated           -0.06662    0.12453 63.00000  -0.535   0.5945
## pHambient                   -0.07462    0.12453 63.00000  -0.599   0.5512
## Temperatureheated:pHambient  0.30557    0.18089 63.00000   1.689   0.0961
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717              
## pHambient   -0.717  0.514       
## Tmprtrhtd:H  0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.10158 0.10158     1    63  0.7492 0.39000  
## pH             0.08185 0.08185     1    63  0.6038 0.44006  
## Temperature:pH 0.38685 0.38685     1    63  2.8534 0.09612 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at tip model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal

## [1] 38  8
shapiro.test(residuals(modTip)) #Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))

# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4545  0.715
##       63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")

Conduct post hoc analysis on calcification at tip.

## POST HOC ANALYSIS 
emmTip <- emmeans(modTip, ~Temperature*pH, adjust = "tukey")  #Estimated marginal means (Least-squares means)
multcomp::cld(emmTip) #create a compact letter display of all pair-wise comparisons
##  Temperature pH        emmean     SE   df lower.CL upper.CL .group
##  ambient     ambient     1.62 0.0868 18.2     1.43     1.80  1    
##  heated      acidified   1.62 0.0868 18.2     1.44     1.81  1    
##  ambient     acidified   1.69 0.0900 18.2     1.50     1.88  1    
##  heated      ambient     1.86 0.0993 16.6     1.65     2.07  1    
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
pwpp(emmTip) #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

ii. Calcification Ratio Analysis at the base

Assemble data for chosen spine bases (proximal end)

### Model 2: calcification in the base of the spines

resultsBase <-
   #create data using only the SEM images of the base of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Base")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results for spine bases with linear mixed model

#LMM for calcification at the tips of spines 
modBase <- 
  resultsBase %>% 
  lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)

summary(modBase) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 98.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5535 -0.6315 -0.1529  0.6918  1.9937 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.1303   0.3610  
##  Residual             0.1597   0.3996  
## Number of obs: 68, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                   2.1440     0.1749 18.8077  12.257 2.06e-10
## Temperatureheated             0.3566     0.2487 19.1707   1.434   0.1677
## pHambient                     0.6524     0.2474 18.8077   2.637   0.0163
## Temperatureheated:pHambient  -0.2224     0.3594 18.9804  -0.619   0.5434
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                   *  
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703              
## pHambient   -0.707  0.497       
## Tmprtrhtd:H  0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Temperature    0.30743 0.30743     1 18.992  1.9251 0.181362   
## pH             1.48873 1.48873     1 18.960  9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114     1 18.980  0.3829 0.543424   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at base model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal

## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))

# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   1.068  0.369
##       64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")

Conduct post hoc analysis on calcification at base

## POST HOC ANALYSIS 
emmBase <- emmeans(modBase, ~Temperature*pH, adjust = "tukey")
  #Estimated marginal means (Least-squares means)
multcomp::cld(emmBase)
##  Temperature pH        emmean    SE   df lower.CL upper.CL .group
##  ambient     acidified   2.14 0.175 18.8     1.78     2.51  1    
##  heated      acidified   2.50 0.177 19.5     2.13     2.87  12   
##  ambient     ambient     2.80 0.175 18.8     2.43     3.16  12   
##  heated      ambient     2.93 0.192 18.8     2.53     3.33   2   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
  #create a compact letter display of all pair-wise comparisons
pwpp(emmBase)

  #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.

iii. Figure 4: Calcification Ratio Interaction Plot at tips and bases

Make Interaction plots for calcification at the base and the tip

Combine both interaction plots at base and tip to be in horizontal line

#Make plots of calcification at tip and base in horizontal line with proper labeling of a and b
#removed legend on the tip, so they can share the same one. But now they are different sizes...

calcificationfig1<-plot_grid(tipplot, baseplot, labels = c("a", "b"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(0.8,1), label_size = 20, label_y=1, align="h");calcificationfig1

#ggsave(filename="Figures/20200611/CalcificationFig1.png", plot=calcificationfig1, dpi=300, width=16, height=6, units="in")

c. Spine Loss

i. Spine Loss Analysis:

Assemble data for dropped spines

## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.

resultsDropped <-
  #create data using only the needed variables.
  results %>% 
  select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>% 
  drop_na(SpineCount)

Analyse dropped spines with linear mixed effect model.

##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <- 
  resultsDropped %>% 
  lm(SpineCount~ Temperature * pH, data=.)

summary(modDropped)
## 
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0000  -6.0833  -0.1667   0.4000  20.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   21.167      3.640   5.814 1.33e-05 ***
## Temperatureheated             -1.167      5.148  -0.227  0.82315    
## pHambient                    -21.000      5.148  -4.079  0.00064 ***
## Temperatureheated:pHambient    1.600      7.461   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
## 
## Response: SpineCount
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## Temperature     1    1.52    1.52  0.0192    0.8914    
## pH              1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH  1    3.66    3.66  0.0460    0.8325    
## Residuals      19 1510.87   79.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for model for dropped spines

##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal

## 1458 1460 
##    2    4
shapiro.test(residuals(modDropped)) #faaaiiilll
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))

# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.8487 0.06483 .
##       19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")

Analysis for dropped spines does not meet assumptions. Attempt a square root transformation.

##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
  # transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
  # run lm model of transformed data

###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal

## 1458 1460 
##    2    4
shapiro.test(residuals(modDropped1)) #fail
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))

# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")

summary(modDropped1)
## 
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0000 -3.0417 -0.0833  0.2000 10.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  10.5833     1.8202   5.814 1.33e-05 ***
## Temperatureheated            -0.5833     2.5742  -0.227  0.82315    
## pHambient                   -10.5000     2.5742  -4.079  0.00064 ***
## Temperatureheated:pHambient   0.8000     3.7304   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
## 
## Response: tdata
##                Sum Sq Df F value    Pr(>F)    
## Temperature      0.23  1  0.0118    0.9146    
## pH             586.44  1 29.4995 3.062e-05 ***
## Temperature:pH   0.91  1  0.0460    0.8325    
## Residuals      377.72 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Transfomraiton (of any usual kind) of dropped spines was not successful, so we use a non parametric Kruskal Wallis Test.

### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  resultsDropped$SpineCount by resultsDropped$Treatment
## Kruskal-Wallis chi-squared = 17.505, df = 3, p-value = 0.0005563

Conduct Dunn Post Hoc Test for dropped spines

## POST HOC Pairwise
dunn <- 
  resultsDropped %>% 
  dunnTest(SpineCount ~ Treatment,
           method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
##   Group Letter MonoLetter
## 1   AMB      a        a  
## 2    OA      b         b 
## 3     T     ac        a c
## 4  T/OA     bc         bc

ii. Figure 5: Spine Loss Interaction Plot

droppedsummary<-results %>% 
  select("Temperature","pH","SpineCount") %>%
  drop_na(SpineCount) %>% 
  group_by(pH, Temperature) %>%
  summarize(mean = mean(SpineCount), se = std.error(SpineCount), N = length(SpineCount))

  dropplot<-ggplot(data=droppedsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
    geom_point(size = 3,  position=position_dodge(0.1)) +
    geom_line(aes(group=pH),  position=position_dodge(0.1)) +
    xlab("Temperature Treatment") +
    ylab("Dropped Spines")+
    geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
    ylim(0,30) +
    theme_bw() + 
    scale_x_discrete(limits=c("ambient", "heated"), 
                     label=c("Ambient", "High"))+
    scale_colour_manual(name = "pH Treatment",
                        values = c("skyblue3", "gray40"),
                        labels = c("Acidified", "Ambient"),
                        guide = guide_legend(reverse = TRUE)) +
    #geom_text(x=1.5, y=11, label="p(pH)<0.001", size=6, color="black") + 
    #geom_text(x=1.5, y=9, label="p(Temperature)=0.91", size=6, color="darkgray") + 
    #geom_text(x=1.5, y=7, label="p(Interaction)=0.83", size=6, color="darkgray") +
    theme_classic()+
    theme(legend.position = "right", 
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14), 
          legend.title = element_text(size=16, face="bold"));dropplot
## Warning: Removed 1 rows containing missing values (geom_errorbar).

#ggsave(filename="Figures/20200611/DroppedSpinesFig.png", plot=dropplot, dpi=300, width=8, height=6, units="in")